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ABSTRACT 

In this paper we report the results of the first ever time-dependent general rela- 
tivistic magnetohydrodynamic simulations of the magnetically dominated monopole 
magnetospheres of black holes. It is found that the numerical solution evolves towards 
a stable steady-state solution which is very close to the corresponding force-free solu- 
tion found by Blandford and Znajek. Contrary to the recent claims, the particle inertia 
does not become dynamically important near the event horizon and the force-free ap- 
proximation provides a proper framework for magnetically dominated magnetospheres 
of black holes. For the first time, our numerical simulations show the development of 
an ultra-relativistic particle wind from a rotating black hole. However, the flow re- 
mains Poynting dominated all the way up to the fast critical point. This suggests that 
the details of the so-called "astrophysical load" , where the electromagnetic energy is 
transferred to particles, may have no effect on the efhciency of the Blandford-Znajek 
mechanism. 
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1 INTRODUCTION 

It is now widely believed that the relativistic jets generated 
in active galactic nuclei, galactic microquasars and, presum- 
ably, during gamma ray bursts are powered by rapidly ro- 
tating black holes. This paradigm is largely based on the 
theoretical results of Blandford and Znajek (1977), who ar- 
gued that, under typical conditions of astrophysical black 
holes, their rotational energy can be efficiently extracted in 
the form of magnetically dominated relativistic wind. 

In fact, the original model of Blandford and Znajek was 
based on the approximation of steady-state force-free degen- 
erate electrodynamics (FFDE, see also Macdonald & Thorne 
1982). In this approximation the dynamical role of magneto- 
spheric plasma is reduced to providing perfect conductivity 
in the space surrounding the black hole, and the particle in- 
ertia is totally neglected. Not all mathematical and physical 
aspects of this electrodynamic model are well understood. 
Similar classical systems, e.g. the Faraday disc, involve the 
poloidal current driven over the surface of a rotating con- 
ductor by a non-electrostatic electromotive force. Moreover, 
the inertia of this conductor is a key factor in the origin 
of this force. Obviously, there is no such conductor in the 
Blandford-Znajek model. Often, this missing element of the 
model is artificially introduced in the form of the so-called 
"membrane" , or the "stretched horizon" , located somewhat 
above the real event horizon (IMacdonald fc Thorne 19821 



iTh orne et al. 198611 . Mathematically, this is consistent with 
the so-called "horizon boundary conditions" employed in the 
Blandford-Znajek model jZnajek 1977| l. 

However, the "Membrane paradigm" cannot always 
provide correct insights into the black hole electrodynamics 
(this would be a great mystery if it did.) Since the membrane 
is very close to the real horizon, the outgoing fast waves 
emitted by the membrane are highly red-shifted. Moreover, 
the membrane cannot at all communicate with the outer 
space by means of Alfven waves because the inner Alfven 
surface is quite far away from the horizon. As the result, 
these waves cannot transport angular momentum and ad- 
just the poloidal electric current of the outgoing wind as 
they do so, for example, in stellar winds of magnetised stars. 
Punsly and Coroniti (1990) used these causality arguments 
to expose the lack of physical clarity in the electrodynamic 
model of Blandford and Znajek (1977) and its representa- 
tion in the membrane paradigm. In fact, they concluded that 
there is no physically meaningful electromotive force in the 
Blandford-Znajek model but only the artificial one which 
has been effectively introduced via Znajek's horizon condi- 
tion. If so, then the Blandford-Znajek model and, hence, its 
steady-state solutions have to be nonphysical. Since these 
FFDE solutions are, in fact, proper mathematical solutions 
their nonphysical nature has to give itself away via instabil- 
ity ( |Punsly fc Coroniti 1990a| [Funsly 20011 . 

In fact, more or less the same criticism was applied to 
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the GRMHD model of black hole magnetospheres developed 
by Phinney (1982) as it also utilized Znajek's horizon con- 
dition. In addition to this criticism, Punsly and Coroniti 
developed a rather attractive alternative GRMHD theory of 
black hole magnetospheres where particle inertia played the 
key role in the extraction of rotational energy of black holes 
jPunsly fc Coroniti 1990E| [Punsly lOOil . While almost ev- 
erywhere in their models the magnetic field was completely 
dominating, in some regions within the black hole ergosphere 
the plasma particles were accelerated to such a high Lorentz 
factor that their inertia could no longer be ignored. This led 
to strong inertially driven electric currents flowing across 
the magnetic field lines, so that the system resembled the 
famous Faraday disc rather closely. 

In spite of the great astrophysical importance of 
the electrodynamic/MHD mechanism of extraction of 
rotational energy of black holes and admirable efforts of 
B. Punsly, the dramatic clash between the Blandford-Znajek 
and Punsly-Coroniti models seemed to remain largely 
unnoticed by the astrophysical community for a whole 
decade. Only very recently, following the developments 
in the theory of force-free electrodynamics (Uchida 1 9971 
IGruzinovl999l IKomissaJOv 2002ll and the impressive 
progress in numerical methods for relativistic astrophysics, 
e.g. llPons et al. 19981 |K omissarov "19991 IKoide et al. 19991 
[Koidc 2003 Koldoba ot al.SOOS' Gammic ct al.2003l 
[Del Zannaet al. 2003 Do Villiers fc Hawley 2003 J , things 
began to change. One of most interesting recent results is 
concerned with the conjecture of Punsly and Coroniti on 
instability of the Blandford-Znajek monopole solution. In 
omissarov 2001b'll this conjecture was tested by means 
of time-dependent general relativistic FFDE simulations 
and was found to be incorrect. Contrary to the appealing 
Punsly-Coroniti causality arguments, this steady-state 
monopole solution of Blandford and Znajek is asymptoti- 
cally stable. There is, however, the question of validity of 
the FFDE approximation which has to be fully investigated 
before reaching any final conclusion. 

There are many examples in physics and astrophysics 
where particular approximations fail to describe certain nat- 
ural phenomena because one or another physical factor ig- 
nored in the approximation becomes important. These ex- 
amples taught us always to check the limits of applicability 
of our theoretical models and never be satisfied by their mere 
self-consistency. The very relevant example can be found in 
the theory of magnetically dominated pulsar winds where 
particle inertia becomes an important factor as a result 
of acceleration by electromagnetic field (Mestel 1999, Sec- 
tion 13.2.3) Other potential examples in the theory of black 
hole magnetospheres are discussed by Punsly (2001, Sec- 
tions 8 and 9). Recently, Punsly (2004) argued that FFDE 
is deficient near the event horizon where "the plasma at- 
tains an infinite inertia in a global sense". Macdonald and 
Thorne ( 1982 ) also anticipated a breakdown of the FFDE 
approximation near the event horizon, though no detailed 
explanations were given. Thus, in the problem of electrody- 
namic/MHD mechanism of extraction of rotational energy 
of black holes we have to verify that the particle inertia can 
indeed be ignored everywhere without oversimplifying this 
problem. This can be done via GRMHD simulations of black 
hole magnetospheres. 

The main goal of this paper is to explore the role of 



particle inertia in the pair-filled monopole magnetospheres 
of black holes by means of such simulations. In particular, we 
need to know whether incorporation of initially small parti- 
cle inertia will eventually lead to strong deviations from the 
corresponding electrodynamic solution both locally, e.g. near 
the event horizon, and globally. By latter we mean changes 
of the global system of electric current, the angular veloc- 
ity of magnetic field lines, and, hence, the efficiency of the 
energy extraction. Since, at present only codes for the per- 
fect relativistic MHD are readily available, the monopole 
magnetic configuration appears to be the most suitable 
one for this purpose. Indeed, this configuration allows a 
dissipation free FFDE solution (Blan dford fc Znajek 1977| 
IKomissarov 200Tbll . Dipolar magnetospheres, like the one 
considered by Koide II2UU3L are less suitable because the 
corresponding electrodynamic solution involves a strong cur- 
rent sheet located within the ergosphere (Komissarov 2002| 
iKomissarov 2004 1 . This indicates that the approximation of 
perfect MHD is also likely to break down in this configura- 
tion. 



2 BASIC EQUATIONS AND NUMERICAL 
METHOD 

In this study we use the Kerr-Schild coordinates, which are 
more suitable for our purpose than the more popular Boyer- 
Lindquist coordinates as they do not introduce a coordinate 
singularity at the event horizon. Not only this allows to in- 
crease efficiency of computer simulations but we may also 
place the inner boundary of our computational domain in- 
side the event horizon. At such a boundary we may confi- 
dently impose usual radiative boundary conditions. Thus, 
the disputed Znajek's horizon conditions are not involved in 
any form. 

Using the standard notation of the 3-1-1 approach, the 

is 



ds 
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where 7^-, is the metric tensor of space. 
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a is the lapse function. 



a = 1/VTTZ, 

and f} is the shift vector, 

= (0,^,0). 

In these equations 



2 2,22, 
K — r + a cos 6 



(1) 



(2) 



(3) 



(4) 



Z = 2r/ti^ 



A = + - 2r. 



and we assume that indexes 1,2,3 correspond to the (p-, r-, 
and ^-coordinates respectively. Notice that we utilize such 
units that G = M = c = 1 and An does not appear in the 
Maxwell equations. 
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The system of perfect general relativistic magnetohy- 
drodynamics (GRMHD) includes the continuity equation, 



dtiy^pu) + di{yf)pu) = 0, 
the energy-momentum equations, 



and the induction equation, 
9t(B') + e'^"a,(£fe) = 0. 



(5) 
(6) 

(7) 



Here gafi is the metric tensor of spacetime, 7 = det(7ij), 
e*-''' is the Levi-Civita pseudo-tensor of space, p is the proper 
mass density of plasma, is its four- velocity vector. The 
total stress-energy- momentum tensor, T'^" , is a sum of the 
stress-energy momentum tensor of matter. 



(8) 



Till' LL V U 

(m) = wu'^u - pg"^ 

where p is the thermodynamic pressure and w is the enthalpy 
per unit volume, and the stress-energy momentum tensor of 
electromagnetic field, 

1 



li)g>^\ (9) 



where F"^ is the Maxwell tensor of the electromagnetic field. 
The electric field, J5, and the magnetic field, i?, are defined 
via 

a 
2' 

and 



(10) 



(11) 



where *F'^'^ is the Faraday tensor of the electromagnetic 
field, which is simply dual to the Maxwell tensor. In the 
limit of ideal MHD 



E 



-vxB or E^^e^jkV^B 



(12) 



where = /u^ is the usual 3- velocity of plasma. All the 
components of vectors and tensors appearing in equations 
Il5l(il71 are the components in the basis of global coordinates. 

Our numerical scheme is a generalization of the scheme 
for the special relativistic MHD described in Komis- 
sarov(1999), where the magnetic field is evolved using the 
method of constraint transport l |Evans fc Hawley 1988| . It 
shares many common features with other existing numeri- 
cal schemes for GRMHD and for this reason we only briefly 
outline its design. 

This is a 2-dimensional scheme with enforced axial sym- 
metry. The set of primitive variables include P, p, the com- 
ponents of magnetic field, B*, and the components of fluid 
four- velocity, as measured in the orthonormal basis of the 
local flducial observer (FIDO), {ii,i^,if,i§} ■ The transfor- 
mation from the coordinate basis to the orthonormal basis 
of FIDO is determined by 

■ - 

* a dt a dr^ 

1 d 



A 



l4><t> d(l> ' 

d jr,p d 
dr d(j) 



(13) 



where A — yT^I^Tee/^. Notice, that p, P, and B'^ are 
deflned at the geometrical centers of the computational cells, 
whereas B*^ and B^ are defined at the geometrical centers 
of the cell interfaces. 

At the beginning of each time step both B'' and B^ 
are found in the cell centers via linear interpolation. The 
cell centered values are then used to set the Riemann prob- 
lems at the cell interfaces using a slope-limited interpolation. 
These problems are then solved using the special relativistic 
Riemann solver described in (IKomissarov 3999J, the solu- 
tion being interpreted as the one observed by the FIDO ini- 
tially located at the interface. In order to find the resolved 
state at the interface, the relative motion of the interface 
and its FIDO has to be taken into account ( Pons ct al.l998|l. 
This resolved state allows to compute the interface fluxes of 
all conserved quantities , first in the FIDO basis and then 
in the coordinate basis via the corresponding transforma- 
tion laws. These fluxes are then used to evolved the volume 
averaged mass, energy, and momentum densities, as well as 
B"^ . The values of in the resolved states at the cell inter- 
faces allow to find the mean values of E^ at the cell edges. 
These are used to evolve the magnetic fiux through the cell 
interfaces and, hence, to update the r- and ^-components 
of the interface magnetic field, first in the coordinate basis 
and then in the FIDO's basis via the corresponding trans- 
formation law. Next, the updated r- and components of 
magnetic field at the cell centers are found via linear inter- 
polation. Finally, the updated conserved variable are trans- 
formed into FIDO's basis and the updated primitive vari- 
ables are found using the same iterative procedure as in 
the special relativistic scheme. The second order accuracy in 
time is achieved via a half time step as described in Komis- 
sarov (1999). 

The scheme was tested using the suit described in 
JKoide et al. 19991 . In all cases, where the test problem was 
based on an existing analytical solution, the agreement be- 
tween numerical and analytical solutions is very good. How- 
ever, the suite also includes the "sub-Keplerian disc" prob- 
lem for which there is no analytical solution. Surprisingly, 
our numerical solution for this problem is dramatically dif- 
ferent from the one of Koide et al.(1999) - instead of "bounc- 
ing of the centrifugal barrier" our disc gets swallowed by 
the hole. Since the specific angular momentum of this disc 
is lower than the one of the last stable orbit, the "bounc- 
ing" seems highly unlikely and we suspect that the re- 
sults by Koide et al.(1999) are incorrect (in fact, this has 
been confirmed in private communication to the author by 
K.Shibata.) In addition, we considered the problem of an 
equilibrium torus around a rotating black hole, both with 
and without an azimuthal magnetic field. The pure gas- 
dynamic solution is described in llAbramowicz et al. 197811 . 
The solution for a magnetized torus does not seem to have 
been described in the literature (a paper is being prepared 
for publication elsewhere.) 
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Figure 1. Solution for a black hole with a = 0.9 at t = 170. Left panel: Thin contour lines show the distribution of the Lorentz factor, 
W. There are 15 contours equally spaced between W = 1 and W = 13.6. Along the equatorial plane W is gradually increasing until it 
reaches a maximum at the cylindrical radius ro ft! 150. The thick lines show the Alfven and fast magnetosonic critical surfaces of the 
outgoing wind. Right panel: The inner region of this solution. The thin lines show the magnetic flux surfaces. The thick lines show, in 
the order of increasing distance from the black hole, the Alfven surface of the ingoing wind, the ergosphere, the wind separation surface, 
and the Alfven critical surface of the outgoing wind. The inner fast surface is too close to the event horizon to be seen in this figure. 



3 NUMERICAL SIMULATIONS 

In these simulations, the computational grid has 100 cells 
in the ^-direction, where it is uniform, and 150 cells in the 
r-direction. The cell size in the r-direction is such that the 
corresponding physical lengths in both directions are equal. 
The radiative outer boundary is located far away from the 
event horizon, rout ~ 230, which ensures that it does not 
effect the solution near the black hole. 

The initial solution describes a purely radial magnetic 
field, 



Bo sin 9 d 
^ 9r' 



(14) 



with Bo = 1 (this implies a non-vanishing azimuthal com- 
ponent in the Boyer-Lindquist basis.) The plasma velocity 
relative to FIDO is set to zero, = 0, whereas its pressure 
and density are set to the same value, of p = p = SVlOO. 
In any conservative scheme for the relativistic magnetohy- 
drodynamics, there is an upper limit on the magnetization 
of plasma above which the hydrodynamical part of the solu- 
tion suffers from large numerical errors (Komissarov 2001a 
IGammie et al.2003t . At this limit, which depends on the 
resolution, the numerical error for the energy density of the 
electromagnetic field becomes comparable with the energy 
density of matter. For this reason we could not set the ini- 
tial density and pressure to the much lower values typical for 
black hole magnetospheres. Moreover, soon after the start of 
the simulations, the numerical solution can be described as 
a pair of particle winds, one ingoing and one outgoing (see 
figure 1^1, and, as the result, the plasma density gradually 



slides down towards the danger zone near, as well as inside, 
of the wind separation interface. To keep the magnetization 
reasonably low, we had to continuously pump new plasma in 
this region. The critical condition we set in these simulations 

was 



«;W^-p = 0.03B^ 



(15) 



where W is the Lorentz factor of the flow as measured 
by FIDO. When the energy density of matter dropped be- 
low 0.03B^ both p and p were artificially increased by 
the same factor. To minimize the effect of the mass in- 
jection on the winds kinematics the velocity of the in- 
jected matter was set to be equal to the local velocity of 
the wind. In fact, new particles must be constantly cre- 
ated in real magnetospheres of black holes but the details 
of this process can be rather different lIBeskin et al. 19911 
IHirotani fc Okamoto IME P hinney 1982| l. 

An additional lower limit was set on the value of the 
thermodynamic pressure, which was not allowed to drop be- 
low O.Olp. In these simulations we used the polytropic equa- 
tion of state with the ratio of specific heats, F = 4/3. 

Figure Q shows the numerical solution for a black hole 
with a = 0.9 a,t t = 170. In the left panel of this figure, 
where the distribution of the Lorentz factor as measured by 
the local FIDO is shown, one can see an almost spherical 
wave front which designates the expanding boundary of the 
outgoing wind. This wind is ultra-relativistic and superfast 
within most of its volume. Moreover, the positions of both 
the fast and the Alfven critical surfaces in the equatorial 
plane do no longer show any noticeable variation at this 
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Figure 2. The angular distribution of _ff^ and fl for a black hole with a = 0.9 at t = 170. Left panel: the diamonds show the MHD 
solution at r = 5 and the crosses show this solution at r = 50. The continuous line show the corresponding FFDE solution at r = 5. 
Right panel: f2; the diamonds show the MHD solution at r = 5 and the crosses show this solution at r = 50. 



point. The wind is magnetically dominated but its magne- 
tization is slowly decreasing with distance. 

While the Lorentz factor of the outgoing wind gradually 
increases with the distance from the black hole, and so does 
the inertia of the accelerated particles, the Lorentz factor 
of the ingoing wind remains lower than = 2 all the way 
down to the inner boundary of the computational domain 
(see also figure Thus, the electromagnetic field does not 
push the plasma of the ingoing wind onto the high- velocity 
orbits and its inertia does not become a key factor in the 
fiow dynamics. On the contrary, new particle have to be 
constantly injected into this wind to keep its inertia above 
the low level given by ea llSI The region of particle injection 
is elongated along the symmetry axis with the major semi- 
axis about 20 and the minor semi-axis about 6. 

The right panel of this figure shows the key sur- 
faces in the inner part of this solution where it has set- 
tled to a steady-state. The comparison with fig. lb in 
l|KomissaJOv 2001bll reveals that the locations of the Alfven 
surfaces are very close to those in the corresponding FFDE 
solution. As expected, the surface separating the ingoing 
wind from the outgoing one is located between the Alfven 
surfaces (Ta kahasln et ai.l990t . 

Figure |5| shows the distribution of the angular velocity 
of magnetic field lines, fi, and the J/^-component of vector 
H defined via 

= f e,,,F^\ (16) 
In steady-state 

VxH^J, (17) 

where J is the electric current density ijKomissarov 2004)1 
and, thus, is a measure of the poloidal current. (In Bland- 
ford & Znajek (1977) i?^ is denoted as Bt-) f2 and are 



very important quantities as they determine the poloidal 
fiuxes of the electromagnetic energy and angular momentum 
(Blandford & Znajek 1977| |Komissaro v 200111 . For example, 
the energy fiux density, Sp, is given by 

Sp = -{H^n)Bp, (18) 

where Bp is the poloidal component of the magnetic field. 
As one can see in figure |21 the distributions of f2 and in 
our magnetically dominated MHD solution arc very close to 
those found in the FFDE solution (Komissarov 2001bL all 
the way up to the fast critical surface of the outgoing wind. 
Thus, the particle inertia does not have much of an effect on 
the process of extraction of energy and angular momentum, 
which remains essentially electromagnetic. 

Figure |3 shows the velocity distribution in the equa- 
torial plane near the black hole. At the point separating 
the infiow from the outflow the angular velocity of plasma, 
= d(f)/dt, equals to the one of the magnetic field lines. 
Closer to the black hole > Q and near the event horizon 
it even exceeds the angular velocity of the black hole, flh- 
(Notice, that it is only the Boyer-Lindquist u*^ which always 
equals to Qh at the event horizon.) 

Figure |3 also shows that the Lorentz factor of the 
ingoing wind remains very low and, thus, there is no 
any reason for the breakdown of the FFDE approxi- 
mation at the event horizon suggested in jPunsly 2004| 
IMacdonald fc Thorne 19821 . 



4 DISCUSSION AND CONCLUSIONS 

In this paper we described the first ever GRMHD simula- 
tions of monopole magnetospheres of rotating black holes. 
Given the central role played by the monopole problem in 



© 0000 RAS, MNRAS 000, 000-000 



6 S.S. Komissarov 



2.0 



1.5 



1.0 



0.5 



0.0 




Figure 3. The Lorentz factor as measured by FIDO, W, the 
radial velocity of plasma, v"^ = dr/dt, the angular velocity of 
plasma, v'f' = d(f>/dt, and the angular velocity of magnetic field 
lines, n, in the equatorial plane at t = 170. 



the development of the general theory of black hole magne- 
tospheres it is not surprising that the results of these simu- 
lations allow to re-examine a number of important issues of 
the theory. 

The main conclusion which follows from these results 
is the validity of the FFDE approximation at least in the 
case of the monopole magnetic configuration considered 
by Blandford and Znajek(1977). Inertia of initially rarefied 
magnetospheric plasma does not grow dynamically impor- 
tant neither in any localized regions near the black hole hori- 
zon nor in a global sense. The system of poloidal electric 
current and the efficiency of energy extraction are basically 
the same as in our GRMHD solution as in the FFDE so- 
lution found earlier lIKomissarov 2001bll . This means that 
in the magnetospheres of black holes there exists an elec- 
tromotive force that has nothing to do with particle iner- 
tia. In a separate paper IKomissarov we show that, in 
great contrast to the Faraday disc or a magnetized stellar 
wind, the poloidal currents in the Blandford-Znajek model 
are driven by the so-called "gravitationally induced" electric 
field, which was first discovered by Wald (1974). Contrary 
to the conclusion reached in ( |Punsly fc Coroniti 1990at , this 
field cannot be screened within the black hole ergosphere by 
any static distribution of electric charge. We also show there 
that Znajek's horizon condition is not a boundary condition 
after all, but a regularity condition imposed at the fast crit- 
ical point of the ingoing wind in the limit of vanishing par- 
ticle inertia. This proves the legitimacy of its utilization in 
the steady-state solutions by Blandford-Znajek(1977) and 
Phinney(1982). 

In addition, the numerical solution shows a number of 
other interesting features that deserve discussing. 

Our GRMHD solution remains very close to the FFDE 



one all the way up to the fast critical surface of the outgoing 
wind and even at the fast critical surface the wind is still 
Poynting fiux dominated. Further away the electromagnetic 
energy may be transferred to the wind particles. However, 
the details of this energy transfer as well as the details of the 
interaction between the wind and its surrounding which is 
responsible for the observational phenomena like superlumi- 
nal jets, radio galaxies etc. cannot effect the wind solution 
in the sub-fast region. This makes us wonder whether the 
key global properties of the Blandford-Znajek mechanism, 
such as its efficiency, are sensitive at all to the nature of the 
so-called "astrophysical load" (^acdonald & Thorn e"l982l 
IThorne et al. 198611 . 

Given the potential importance of this finding let us 
show that a RMHD flow may indeed stay Poynting flux dom- 
inated at the fast critical point. For this purpose we consider 
a one dimensional flow in flat spacetime. In the limit of cold 
MHD, the ratio kinetic energy flux to the Poynting fiux is 

K = pW'/Bl 

where Bt is the tangential component of magnetic field. At 
the fast point the fiow velocity 

where 6^ = + (Bt/W)'^ , B„ being the normal component 
of the magnetic field. From these two equations one finds 
that at the fast point 



1 _B| 



(19) 



Thus, the fiow is Poynting dominated provided the fast 
speed is ultra-relativistic and the magnetic field is predom- 
inantly tangential. 

These are the first GRMHD simulations where an 
ultrarelativistic outfiow is produced by a black hole. In 
all previous simulations of this sort (]Roidc ct al. 19991 
IKoide et a l?2000 Koidc 200 ^ Ikomissarov 2001al 

IGammie^t al.2003 ,De Villiers &i Hawley 2003| l no such 
high velocity outfiows were reported, which was a bit 
worrying given the well known results of various astro- 
physical observations. In more realistic simulations of the 
future, which will include both the accretion disc and its 
magnetised corona, it may be possible to obtain not only 
ultra-relativistic but also well collimated outflows. 

Since the fast and the Alfven surfaces of the outgoing 
wind are relatively far away from the event horizon, one 
may expect the outgoing wind to be more or less accurately 
described by the fiat space monopole solution. According 
to lIBeskin 1997t the ratio of the fast suface radius, rj, and 
the Alfven surface radius, ra, in the equatorial plane of this 
solution is 



Tf/ra 



.1/3 



whereas the Lorentz factor at the fast surface 



Wf 



.1/3 



where a is the magnetization papameter. In our simulations 
(7^/=* « 4.08, Tf/ra = 6.60, and Wf = 4.35. Thus, the nu- 
merical and the analytical results agree quite well. 

The presented results also suggest a way of dealing with 
such more realistic problems of computational astrophysics 
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of black holes which involve both magnetically dominated ul- 
trarelativistic components and particle dominated slow com- 
ponents. Namely, one may use the same GRMHD approxi- 
mation to model all components of the system but impose 
an upper limit on the magnetization of the ultrarelativistic 
component, via appropriate floor values for both the pres- 
sure and density of matter. Although, there is nothing par- 
ticularly new in such an approach and utilization of similar 
lower/upper limits is widely spread in computational prac- 
tice, the accurate representation of magnetically dominated 
component in our simulations is reassuring. 
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